home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / curvefit.pro < prev    next >
Text File  |  1997-07-08  |  12KB  |  267 lines
  1. ; $Id: curvefit.pro,v 1.13 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1982-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;
  6. function curvefit, x, y, weights, a, sigma, Function_Name = Function_Name, $
  7.                         itmax=itmax, iter=iter, tol=tol, chi2=chi2, $
  8.                         noderivative=noderivative,chisq=chisq
  9. ; Copyright (c) 1988-1995, Research Systems, Inc.  All rights reserved.
  10. ;       Unauthorized reproduction prohibited.
  11. ;+
  12. ; NAME:
  13. ;       CURVEFIT
  14. ;
  15. ; PURPOSE:
  16. ;       Non-linear least squares fit to a function of an arbitrary 
  17. ;       number of parameters.  The function may be any non-linear 
  18. ;       function.  If available, partial derivatives can be calculated by 
  19. ;       the user function, else this routine will estimate partial derivatives
  20. ;       with a forward difference approximation.
  21. ;
  22. ; CATEGORY:
  23. ;       E2 - Curve and Surface Fitting.
  24. ;
  25. ; CALLING SEQUENCE:
  26. ;       Result = CURVEFIT(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
  27. ;                         ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)
  28. ;
  29. ; INPUTS:
  30. ;       X:  A row vector of independent variables.  This routine does
  31. ;        not manipulate or use values in X, it simply passes X
  32. ;        to the user-written function.
  33. ;
  34. ;       Y:  A row vector containing the dependent variable.
  35. ;
  36. ;  Weights:  A row vector of weights, the same length as Y.
  37. ;            For no weighting,
  38. ;                 Weights(i) = 1.0.
  39. ;            For instrumental (Gaussian) weighting,
  40. ;                 Weights(i)=1.0/sigma(i)^2
  41. ;            For statistical (Poisson)  weighting,
  42. ;                 Weights(i) = 1.0/y(i), etc.
  43. ;
  44. ;       A:  A vector, with as many elements as the number of terms, that 
  45. ;           contains the initial estimate for each parameter.  If A is double-
  46. ;           precision, calculations are performed in double precision, 
  47. ;           otherwise they are performed in single precision. Fitted parameters
  48. ;           are returned in A.
  49. ;
  50. ; KEYWORDS:
  51. ;       FUNCTION_NAME:  The name of the function (actually, a procedure) to 
  52. ;       fit.  If omitted, "FUNCT" is used. The procedure must be written as
  53. ;       described under RESTRICTIONS, below.
  54. ;
  55. ;       ITMAX:  Maximum number of iterations. Default = 20.
  56. ;       ITER:   The actual number of iterations which were performed
  57. ;       TOL:    The convergence tolerance. The routine returns when the
  58. ;               relative decrease in chi-squared is less than TOL in an 
  59. ;               interation. Default = 1.e-3.
  60. ;       CHI2:   The value of chi-squared on exit (obselete)
  61. ;     
  62. ;       CHISQ:   The value of reduced chi-squared on exit
  63. ;       NODERIVATIVE:   If this keyword is set then the user procedure will not
  64. ;               be requested to provide partial derivatives. The partial
  65. ;               derivatives will be estimated in CURVEFIT using forward
  66. ;               differences. If analytical derivatives are available they
  67. ;               should always be used.
  68. ;
  69. ; OUTPUTS:
  70. ;       Returns a vector of calculated values.
  71. ;       A:  A vector of parameters containing fit.
  72. ;
  73. ; OPTIONAL OUTPUT PARAMETERS:
  74. ;       Sigmaa:  A vector of standard deviations for the parameters in A.
  75. ;
  76. ; COMMON BLOCKS:
  77. ;       NONE.
  78. ;
  79. ; SIDE EFFECTS:
  80. ;       None.
  81. ;
  82. ; RESTRICTIONS:
  83. ;       The function to be fit must be defined and called FUNCT,
  84. ;       unless the FUNCTION_NAME keyword is supplied.  This function,
  85. ;       (actually written as a procedure) must accept values of
  86. ;       X (the independent variable), and A (the fitted function's
  87. ;       parameter values), and return F (the function's value at
  88. ;       X), and PDER (a 2D array of partial derivatives).
  89. ;       For an example, see FUNCT in the IDL User's Libaray.
  90. ;       A call to FUNCT is entered as:
  91. ;       FUNCT, X, A, F, PDER
  92. ; where:
  93. ;       X = Variable passed into CURVEFIT.  It is the job of the user-written
  94. ;        function to interpret this variable.
  95. ;       A = Vector of NTERMS function parameters, input.
  96. ;       F = Vector of NPOINT values of function, y(i) = funct(x), output.
  97. ;       PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
  98. ;               PDER(I,J) = DErivative of function at ith point with
  99. ;               respect to jth parameter.  Optional output parameter.
  100. ;               PDER should not be calculated if the parameter is not
  101. ;               supplied in call. If the /NODERIVATIVE keyword is set in the
  102. ;               call to CURVEFIT then the user routine will never need to
  103. ;               calculate PDER.
  104. ;
  105. ; PROCEDURE:
  106. ;       Copied from "CURFIT", least squares fit to a non-linear
  107. ;       function, pages 237-239, Bevington, Data Reduction and Error
  108. ;       Analysis for the Physical Sciences.  This is adapted from:
  109. ;    Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear
  110. ;    Parameters", J. Soc. Ind. Appl. Math., Vol 11, no. 2, pp. 431-441,
  111. ;    June, 1963.
  112. ;
  113. ;       "This method is the Gradient-expansion algorithm which
  114. ;       combines the best features of the gradient search with
  115. ;       the method of linearizing the fitting function."
  116. ;
  117. ;       Iterations are performed until the chi square changes by
  118. ;       only TOL or until ITMAX iterations have been performed.
  119. ;
  120. ;       The initial guess of the parameter values should be
  121. ;       as close to the actual values as possible or the solution
  122. ;       may not converge.
  123. ;
  124. ; EXAMPLE:  Fit a function of the form f(x) = a * exp(b*x) + c to
  125. ;    sample pairs contained in x and y.
  126. ;    In this example, a=a(0), b=a(1) and c=a(2).
  127. ;    The partials are easily computed symbolicaly:
  128. ;        df/da = exp(b*x), df/db = a * x * exp(b*x), and df/dc = 1.0
  129. ;
  130. ;        Here is the user-written procedure to return F(x) and
  131. ;        the partials, given x:
  132. ;       pro gfunct, x, a, f, pder    ; Function + partials
  133. ;      bx = exp(a(1) * x)
  134. ;         f= a(0) * bx + a(2)        ;Evaluate the function
  135. ;         if N_PARAMS() ge 4 then $    ;Return partials?
  136. ;        pder= [[bx], [a(0) * x * bx], [replicate(1.0, N_ELEMENTS(f))]]
  137. ;       end
  138. ;
  139. ;         x=findgen(10)            ;Define indep & dep variables.
  140. ;         y=[12.0, 11.0,10.2,9.4,8.7,8.1,7.5,6.9,6.5,6.1]
  141. ;         Weights=1.0/y            ;Weights
  142. ;         a=[10.0,-0.1,2.0]        ;Initial guess
  143. ;         yfit=curvefit(x,y,Weights,a,sigma,function_name='gfunct')
  144. ;      print, 'Function parameters: ', a
  145. ;         print, yfit
  146. ;       end
  147. ;
  148. ; MODIFICATION HISTORY:
  149. ;       Written, DMS, RSI, September, 1982.
  150. ;       Does not iterate if the first guess is good.  DMS, Oct, 1990.
  151. ;       Added CALL_PROCEDURE to make the function's name a parameter.
  152. ;              (Nov 1990)
  153. ;       12/14/92 - modified to reflect the changes in the 1991
  154. ;            edition of Bevington (eq. II-27) (jiy-suggested by CreaSo)
  155. ;       Mark Rivers, U of Chicago, Feb. 12, 1995
  156. ;           - Added following keywords: ITMAX, ITER, TOL, CHI2, NODERIVATIVE
  157. ;             These make the routine much more generally useful.
  158. ;           - Removed Oct. 1990 modification so the routine does one iteration
  159. ;             even if first guess is good. Required to get meaningful output
  160. ;             for errors. 
  161. ;           - Added forward difference derivative calculations required for 
  162. ;             NODERIVATIVE keyword.
  163. ;           - Fixed a bug: PDER was passed to user's procedure on first call, 
  164. ;             but was not defined. Thus, user's procedure might not calculate
  165. ;             it, but the result was then used.
  166. ;
  167. ;    Steve Penton, RSI, June 1996.
  168. ;            - Changed SIGMAA to SIGMA to be consistant with other fitting 
  169. ;              routines.
  170. ;            - Changed CHI2 to CHISQ to be consistant with other fitting 
  171. ;              routines.
  172. ;            - Changed W to Weights to be consistant with other fitting 
  173. ;              routines.
  174. ;         _ Updated docs regarding weighing.
  175. ;           
  176. ;-
  177.        on_error,2              ;Return to caller if error
  178.  
  179.        ;Name of function to fit
  180.        if n_elements(function_name) le 0 then function_name = "FUNCT"
  181.        if n_elements(tol) eq 0 then tol = 1.e-3        ;Convergence tolerance
  182.        if n_elements(itmax) eq 0 then itmax = 20    ;Maximum # iterations
  183.     type = size(a)
  184.     type = type[type[0]+1]
  185.     double = type eq 5
  186.     if (type ne 4) and (type ne 5) then a = float(a)  ;Make params floating
  187.  
  188.        ; If we will be estimating partial derivatives then compute machine
  189.        ; precision
  190.        if keyword_set(NODERIVATIVE) then begin
  191.           res = machar(DOUBLE=double)
  192.           eps = sqrt(res.eps)
  193.        endif
  194.  
  195.        nterms = n_elements(a)   ; # of parameters
  196.        nfree = n_elements(y) - nterms ; Degrees of freedom
  197.        if nfree le 0 then message, 'Curvefit - not enough data points.'
  198.        flambda = 0.001          ;Initial lambda
  199.        diag = lindgen(nterms)*(nterms+1) ; Subscripts of diagonal elements
  200.  
  201. ;      Define the partial derivative array
  202.        if double then pder = dblarr(n_elements(y), nterms) $
  203.     else pder = fltarr(n_elements(y), nterms)
  204. ;
  205.        for iter = 1, itmax do begin   ; Iteration loop
  206.  
  207. ;         Evaluate alpha and beta matricies.
  208.           if keyword_set(NODERIVATIVE) then begin
  209. ;            Evaluate function and estimate partial derivatives
  210.              call_procedure, Function_name, x, a, yfit
  211.              for term=0, nterms-1 do begin
  212.                 p = a       ; Copy current parameters
  213.                 ; Increment size for forward difference derivative
  214.                 inc = eps * abs(p[term])    
  215.                 if (inc eq 0.) then inc = eps
  216.                 p[term] = p[term] + inc
  217.                 call_procedure, function_name, x, p, yfit1
  218.                 pder[0,term] = (yfit1-yfit)/inc
  219.              endfor
  220.           endif else begin
  221.              ; The user's procedure will return partial derivatives
  222.              call_procedure, function_name, x, a, yfit, pder 
  223.           endelse
  224.  
  225.           if nterms eq 1 then pder = reform(pder, n_elements(y), 1)
  226.  
  227.           beta = (y-yfit)*Weights # pder
  228.           alpha = transpose(pder) # (Weights # (fltarr(nterms)+1)*pder)
  229.           chisq1 = total(Weights*(y-yfit)^2)/nfree ; Present chi squared.
  230.  
  231.                 ; If a good fit, no need to iterate
  232.       all_done = chisq1 lt total(abs(y))/1e7/NFREE
  233. ;
  234. ;         Invert modified curvature matrix to find new parameters.
  235.  
  236.           repeat begin
  237. ;             c = sqrt(alpha(diag) # alpha(diag))   ;Old way
  238.          c = sqrt(alpha[diag])    ;Minimize chances of over/underflow.
  239.          c = c # c
  240.              array = alpha / c
  241.          array[diag] = array[diag]*(1.+flambda) 
  242.              if (size(array))[(size(array))[0] + 2] eq 1 then $
  243.              array = (1.0 / array) else $ 
  244.              array = invert(array)
  245.              b = a + array/c # transpose(beta) ; New params
  246.              call_procedure, function_name, x, b, yfit  ; Evaluate function
  247.              chisq = total(Weights*(y-yfit)^2)/nfree         ; New chisq Reduced
  248.          if all_done then goto, done
  249.              flambda = flambda*10.                      ; Assume fit got worse
  250.           endrep until chisq le chisq1
  251. ;
  252.           flambda = flambda/100.  ; Decrease flambda by factor of 10
  253.           a=b                     ; Save new parameter estimate.
  254.           if ((chisq1-chisq)/chisq1) le tol then goto,done  ; Finished?
  255.        endfor                        ;iteration loop
  256. ;
  257.        message, 'Failed to converge', /INFORMATIONAL
  258. ;
  259. done:  sigma = sqrt(array[diag]/alpha[diag]) ; Return sigma's
  260.        ;
  261.        ; chi2 is an obselete keyword
  262.        ; chisq is more conistant
  263.        ;
  264.        chi2 = chisq                          ; Return chi-squared
  265.        return,yfit              ;return result
  266. END
  267.